#include "viewport.h"
#include "../math/functions.h"
#include "../math/functionsse.h"
#include <mmintrin.h>

#if WALD
// ================================================================================================
// Get the projected verts
// ================================================================================================
Vector4 Viewport::GetProjectedVerts(const Vector4& vec) {
    float absx = fabs(vec.x);
    float absy = fabs(vec.y);
    float absz = fabs(vec.z);

    uint k;
    Vector4 kvec;
    if (absx > absy) {
        if (absx > absz) {
            kvec = Vector4(1.0f, 0.0f, 0.0f); /* X */ 
            k = 0;
        }
        else {
            kvec = Vector4(0.0f, 0.0f, 1.0f); /* Z */
            k = 2;
        }
    }
    else {
        if (absy > absz)  {
            kvec = Vector4(0.0f, 1.0f, 0.0f); /* Y */ 
            k = 1;
        }
        else {
            kvec = Vector4(0.0f, 0.0f, 1.0f); /* Z */
            k = 2;
        }
    }

    kvec[3] = float(k);

    return kvec;
}

// -- precomputed possible projected dimensions
static uint modulos[] = {0,1,2,0,1};
// ================================================================================================
// Intersection test for a triangle surface
// ================================================================================================
int Viewport::SurfaceIntersection(Ray& ray, float& tmax, const TriangleAcc& triaccel) {
#define ku modulos[triaccel.k+1]
#define kv modulos[triaccel.k+2]

    // -- most significant bit
    BitComparison beta;
    BitComparison gamma;
    BitComparison theta;

    if(DotProduct(ray.direction, triaccel.normal) > 0.0f) 
        return 0;

    const float denom = 1.0f/(ray.direction[triaccel.k] + triaccel.n_u * ray.direction[ku] + triaccel.n_v *
        ray.direction[kv]);
    const float tplane = (triaccel.n_d - ray.origin[triaccel.k] - triaccel.n_u * ray.origin[ku] - triaccel.n_v * 
        ray.origin[kv]) * denom;

    if(!(tplane < tmax && tplane > 0.0f))
        return 0;

    // -- its a hit store the current scalar value
    const float hu = ray.origin[ku] + tplane * ray.direction[ku];
    const float hv = ray.origin[kv] + tplane * ray.direction[kv];

    beta.realpart = hu * triaccel.b_nu + hv * triaccel.b_nv + triaccel.b_d;
    if(beta.intpart & eMostSig)
        return 0;

    gamma.realpart = hu * triaccel.c_nu + hv * triaccel.c_nv + triaccel.c_d;
    if(gamma.intpart & eMostSig)
        return 0;

    theta.realpart = beta.realpart + gamma.realpart - 1.0f;
    if(theta.intpart & eMostSig) {
        tmax = tplane;
        return 1;
    }

    return 0;
}

#endif

#if HAVEL
// ================================================================================================
// Intersection of triangle and ray, except using the havel and herout method
// ================================================================================================
Vector4 Viewport::GetProjectedVerts(const Vector4& vec){return vec;}
int Viewport::SurfaceIntersection(Ray& ray, float& tmax, const TriangleAcc& triaccel) {
    #ifdef CPP
        Vector4 normal(triaccel.nx, triaccel.ny, triaccel.nz);
        Vector4 norm1(triaccel.ux, triaccel.uy, triaccel.uz);
        Vector4 norm2(triaccel.vx, triaccel.vy, triaccel.vz);

        float det = DotProduct(ray.direction, normal);
        float dett = triaccel.nd - DotProduct(ray.origin, normal);

        float val = dett - ray.direction[3] * det;

        if(*reinterpret_cast<int*>(&val) ^ *reinterpret_cast<int*>(&dett))
            return 0;

        Vector4 detp = (det * ray.origin) + (dett * ray.direction);
        float detu = DotProduct(detp, norm1) + (det * triaccel.ud);
        float u = detu/det;
        if(u < 0.0f || u > 1.0f)
            return 0;
        float detv = DotProduct(detp, norm2) + (det * triaccel.vd);
        float v = detv/det;
        if(v < 0.0f || v > 1.0f)
            return 0;
        if((u + v) > 1.0f )
            return 0;

        ray.origin[3] = dett/det;

        return 1;


    #else

        const float _declspec(align(16)) int_coeff_arr[4] = {-1, -1, -1, 1};
        const float _declspec(align(16)) int_coeff_arr1[4] = {1, 0, 0, 0};
        const __m128 int_coef = _mm_load_ps(&(int_coeff_arr[0]));
        const __m128 int_coef1 = _mm_load_ps(&int_coeff_arr1[0]);

        const __m128 o = _mm_load_ps(&ray.origin[0]);
        const __m128 d = _mm_load_ps(&ray.direction[0]);
        const __m128 n = _mm_load_ps(&triaccel.nx);

        if(!_mm_movemask_ps(DotProductSSE(d, n))) 
            return 0;

        const __m128 det = DotProductSSE(n, d);
        const __m128 dett = DotProductSSE( _mm_mul_ps(int_coef, n), o);
        const __m128 oldt = _mm_load_ss(&tmax);

        if((_mm_movemask_ps(_mm_xor_ps(dett, _mm_sub_ss(_mm_mul_ss(oldt, det), dett)))&1) == 0) {
            const __m128 detp = _mm_add_ps(_mm_mul_ps(o, det), _mm_mul_ps(dett, d));
            const __m128 detu = DotProductSSE(detp, _mm_load_ps(&triaccel.ux));

            if((_mm_movemask_ps(_mm_xor_ps(detu, _mm_sub_ss(det, detu)))&1) == 0) {
                const __m128 detv = DotProductSSE(detp, _mm_load_ps(&triaccel.vx));

                if((_mm_movemask_ps(_mm_xor_ps(detv, _mm_sub_ss(det, _mm_add_ss(detu,detv))))&1) == 0) {
                    const __m128 inv_det = _mm_div_ps(int_coef1, det);
                    _mm_store_ss(&tmax, _mm_mul_ss(dett, inv_det));
                    return 1;
                }
            }
        }

        return 0;

    #endif
};

// -- cache a mask lookup
static const float _declspec(align(16)) validate[16][4] = {1.0f, 1.0f, 1.0f, 1.0f,
                                                           0.0f, 1.0f, 1.0f, 1.0f,
                                                           1.0f, 0.0f, 1.0f, 1.0f,
                                                           0.0f, 0.0f, 1.0f, 1.0f,
                                                           1.0f, 1.0f, 0.0f, 1.0f,
                                                           0.0f, 1.0f, 0.0f, 1.0f,
                                                           1.0f, 0.0f, 0.0f, 1.0f,
                                                           0.0f, 0.0f, 0.0f, 1.0f,
                                                           1.0f, 1.0f, 1.0f, 0.0f,
                                                           0.0f, 1.0f, 1.0f, 0.0f,
                                                           1.0f, 0.0f, 1.0f, 0.0f,
                                                           0.0f, 0.0f, 1.0f, 0.0f,
                                                           1.0f, 1.0f, 0.0f, 0.0f,
                                                           0.0f, 1.0f, 0.0f, 0.0f,
                                                           1.0f, 0.0f, 0.0f, 0.0f,
                                                           0.0f, 0.0f, 0.0f, 0.0f};

static const float _declspec(align(16)) intcoeffarr[4][4] = {-1.0f, -1.0f, -1.0f, -1.0f,
                                                             -1.0f, -1.0f, -1.0f, -1.0f,
                                                             -1.0f, -1.0f, -1.0f, -1.0f,
                                                              1.0f,  1.0f,  1.0f,  1.0f};

static const float _declspec(align(16)) safecheck[4] = {0.0f, 0.0f, 0.0f, 0.0f};
static const float _declspec(align(16)) delta[4] = {0.000001f, 0.000001f, 0.000001f, 0.000001f};

// ================================================================================================
// Surface Intersection with havel and herout's method using data parallelism
// ================================================================================================
void Viewport::SurfaceIntersection(RaySSE& rays, float tmax[], const TriangleAccDP& triaccel) {
    uint uvmask;
    __m128 intcoef[4];
    __m128 oldt;
    __m128 dett0[4];
    __m128 detp0[4];  
    __m128 detp1[4];  
    __m128 detp[4];  
    __m128 detu;
    __m128 detv;
    __m128 o[4];
    __m128 d[4];
    __m128 n[4];
    __m128 u[4];
    __m128 v[4];
    __m128 safe = _mm_load_ps(safecheck);
    __m128 deltas = _mm_load_ps(delta);
    intcoef[0] = _mm_load_ps(intcoeffarr[0]);
    intcoef[1] = _mm_load_ps(intcoeffarr[1]);
    intcoef[2] = _mm_load_ps(intcoeffarr[2]);
    intcoef[3] = _mm_load_ps(intcoeffarr[3]);

    // -- 4x4 of ray origin
    o[0] = _mm_load_ps(rays.ox);
    o[1] = _mm_load_ps(rays.oy);
    o[2] = _mm_load_ps(rays.oz);
    o[3] = _mm_load_ps(rays.ow);

    // -- 4x4 of direction origin
    d[0] = _mm_load_ps(rays.dx);
    d[1] = _mm_load_ps(rays.dy);
    d[2] = _mm_load_ps(rays.dz);
    d[3] = _mm_load_ps(rays.dw);

    // -- 4x4 of plane n
    n[0] = _mm_load_ps(triaccel.nx);
    n[1] = _mm_load_ps(triaccel.ny);
    n[2] = _mm_load_ps(triaccel.nz);
    n[3] = _mm_load_ps(triaccel.nd);

    // -- 4x4 of plane n1
    u[0] = _mm_load_ps(triaccel.ux);
    u[1] = _mm_load_ps(triaccel.uy);
    u[2] = _mm_load_ps(triaccel.uz);
    u[3] = _mm_load_ps(triaccel.ud);

    // -- 4x4 of plane n2
    v[0] = _mm_load_ps(triaccel.vx);
    v[1] = _mm_load_ps(triaccel.vy);
    v[2] = _mm_load_ps(triaccel.vz);
    v[3] = _mm_load_ps(triaccel.vd);

    // -- load the 4 t values into the SSE registers
    oldt = _mm_load_ps(tmax);

    const __m128 det = DotProductSSE4x4(n, d);
    if(_mm_movemask_ps(_mm_cmpge_ps(det, safe)))
        return;

    const __m128 recdet = _mm_rcp_ps(det);

    MulSSE4x4(intcoef, n, dett0);
    const __m128 dett = DotProductSSE4x4(dett0, o);
  

    MulSSE1x4(det, o, detp0);
    MulSSE1x4(dett, d, detp1);
    AddSSE4x4(detp0, detp1, detp);

    detu = DotProductSSE4x4(detp, u);
    detv = DotProductSSE4x4(detp, v);

    // -- check to see if the t is not greater than the max
    uvmask = _mm_movemask_ps(_mm_xor_ps(dett, _mm_sub_ps(_mm_mul_ps(oldt, det), dett)));

    if(uvmask == 15)
        return;

    // -- check to see if u is in the allowable range
    uvmask |=  _mm_movemask_ps(_mm_xor_ps(detu, _mm_sub_ps(det, detu)));

    if(uvmask == 15)
        return;

    // -- check to see if u and v are in the allowable range
    // -- however if we are straddling the edge then allow it
    __m128 rem =  _mm_sub_ps(det, _mm_add_ps(detu, detv));
    uvmask |= (_mm_movemask_ps(_mm_xor_ps(detv, rem)) & 
               _mm_movemask_ps(_mm_cmpnle_ps(FAbsSSE(rem), deltas)));
    if(uvmask == 15)
        return;

    // -- lookup the mask to hide the invalid t values
    __m128 nonvalid = _mm_load_ps(validate[uvmask]);
    _mm_store_ps(rays.t, _mm_mul_ps(_mm_mul_ps(dett, recdet), nonvalid));
}

#endif

#define DP
#ifndef DP
// ================================================================================================
// Write to the whole buffer
// ================================================================================================
void Viewport::Render(MeshObject& mobject, Rect2D& rect) {
    DynamicArray<TriangleAcc> meshacc = mobject.GetTriAccels();
    uint numfaces = mobject.GetNumFaces();
    Ray ray;
    ray.SetOrigin(0.0f, 0.0f, 0.0f);
    
    uint xstart = rect.topleft.x;
    uint ystart = rect.topleft.y;
    uint xend = rect.bottomright.x;
    uint yend = rect.bottomright.y;

    // -- transform the screenspace coordinates to cartesian coordinates
    // -- then compute the steps
    float xstep = (1.0f/float(width)) * (-2.0f);
    float ystep = (1.0f/float(height)) * (-2.0f);
    float xcoordstart = (xstart)/float(width) * (-2.0f) + (1.0f);
    float ycoordstart = (ystart)/float(height) * (-2.0f) + (1.0f);
    unsigned int index = xstart + (ystart * width);
    float aspect = float(width)/float(height);

    // -- iterate the scene
    float ycoord = ycoordstart;
    float tmax = 0.0f;
    for(unsigned int k = ystart; k < yend; ++k) {
        float xcoord = xcoordstart;
        for(unsigned int j = xstart; j < xend; ++j) {
            tmax = 6.0f;
            for(unsigned int i = 0; i < numfaces; ++i) {
                ray.SetDirection(xcoord * aspect, ycoord, depth, 0.0f);

                // check to see if this coordinate needs to be lighted up
                if(SurfaceIntersection(ray, tmax, meshacc[i])) {
                    pixelbuffer[index].r = 1.0f;//(*mobject)[i].GetColor().r;
                    pixelbuffer[index].g = 0.0f;//(*mobject)[i].GetColor().g;
                    pixelbuffer[index].b = 5.0f;//(*mobject)[i].GetColor().b;
                }
            }
            xcoord += xstep;
            ++index;
        }
        ycoord += ystep;
    }
}

#else
// ================================================================================================
// Write to the whole buffer
// ================================================================================================
void Viewport::Render(MeshObject& mobject, Rect2D& rect) {
    DynamicArray<TriangleAccDP> meshacc = mobject.GetTriAccelsDP();
    Ray ray;
    ray.SetOrigin(0.0f, 0.0f, 0.0f);
    uint numfaces = mobject.GetNumFaces();

    RaySSE rays;
    _declspec(align(16)) float tmax[4];
    for(uint i = 0; i < 4; ++i) {
        rays.ox[i] = 0.0f;
        rays.oy[i] = 0.0f;
        rays.oz[i] = 0.0f;
        rays.ow[i] = 1.0f;
        rays.dz[i] = depth;
        rays.dw[i] = 0.0f;
        rays.t[i] = 0.0f;
        tmax[i] = 6.0f;
    }

    uint xstart = rect.topleft.x;
    uint ystart = rect.topleft.y;
    uint xend = rect.bottomright.x;
    uint yend = rect.bottomright.y;

    // -- transform the screenspace coordinates to cartesian coordinates
    // -- then compute the steps
    float xstep = (1.0f/float(width)) * (-2.0f);
    float ystep = (1.0f/float(height)) * (-2.0f);
    float xcoordstart = (xstart)/float(width) * (-2.0f) + (1.0f);
    float xcoord = xcoordstart;
    float ycoordstart = (ystart)/float(height) * (-2.0f) + (1.0f);
    uint index = xstart + (ystart * width);
    float aspect = float(width)/float(height);

    // -- iterate the scene
    float ycoord = ycoordstart;
    float dif;
    float t;
    Vector4 norm;
    for(uint k = ystart; k < yend; ++k) {
        rays.dy[0] = ycoord;
        rays.dy[1] = ycoord;
        rays.dy[2] = ycoord;
        rays.dy[3] = ycoord;

        xcoord = xcoordstart;
        for(uint j = xstart; j < xend; j += 4) {
            rays.dx[0] = xcoord;// * aspect;
            rays.dx[1] = (xcoord + xstep * 1);// * aspect;
            rays.dx[2] = (xcoord + xstep * 2);// * aspect;
            rays.dx[3] = (xcoord + xstep * 3);// * aspect;

            tmax[0] = 6.0f;
            tmax[1] = 6.0f;
            tmax[2] = 6.0f;
            tmax[3] = 6.0f;

            uint i = 0;
            for(; i < numfaces; ++i) {
                rays.t[0] = 0.0f;
                rays.t[1] = 0.0f;
                rays.t[2] = 0.0f;
                rays.t[3] = 0.0f;

                // check to see if this coordinate needs to be lighted up
                SurfaceIntersection(rays, tmax, meshacc[i]);
                norm = Normalize(meshacc[i].nx[0], meshacc[i].ny[0], meshacc[i].nz[0]);
                if(rays.t[0]) {
                    tmax[0] = rays.t[0];
                    dif = fabs(DotProduct(norm, rays.dx[0], rays.dy[0], rays.dz[0]));
                    t = 1.0f/(rays.t[0] * rays.t[0] * rays.t[0]);
            
                    pixelbuffer[index].r = 1.0f * dif * t;//(*mobject)[i].GetColor().r;
                    pixelbuffer[index].g = 0.0f * dif * t;//(*mobject)[i].GetColor().g;
                    pixelbuffer[index].b = 0.0f * dif * t;//(*mobject)[i].GetColor().b;
                }

                if(rays.t[1]) {
                    tmax[1] = rays.t[1];
                    dif = fabs(DotProduct(norm, rays.dx[1], rays.dy[1], rays.dz[1]));
                    t = 1.0f/(rays.t[1] * rays.t[1] * rays.t[1]);
            
                    pixelbuffer[index + 1].r = 1.0f * dif * t;//(*mobject)[i].GetColor().r;
                    pixelbuffer[index + 1].g = 0.0f * dif * t;//(*mobject)[i].GetColor().g;
                    pixelbuffer[index + 1].b = 0.0f * dif * t;//(*mobject)[i].GetColor().b;
                }

                if(rays.t[2]) {
                    tmax[2] = rays.t[2];
                    dif = fabs(DotProduct(norm, rays.dx[2], rays.dy[2], rays.dz[2]));
                    t = 1.0f/(rays.t[2] * rays.t[2] * rays.t[2]);

                    pixelbuffer[index + 2].r = 1.0f * dif * t;//(*mobject)[i].GetColor().r;
                    pixelbuffer[index + 2].g = 0.0f * dif * t;//(*mobject)[i].GetColor().g;
                    pixelbuffer[index + 2].b = 0.0f * dif * t;//(*mobject)[i].GetColor().b;
                }

                if(rays.t[3]) {
                    tmax[3] = rays.t[3];
                    dif = fabs(DotProduct(norm, rays.dx[3], rays.dy[3], rays.dz[3]));
                    t = 1.0f/(rays.t[3] * rays.t[3] * rays.t[3]);

                    pixelbuffer[index + 3].r = 1.0f * dif * t;//(*mobject)[i].GetColor().r;
                    pixelbuffer[index + 3].g = 0.0f * dif * t;//(*mobject)[i].GetColor().g;
                    pixelbuffer[index + 3].b = 0.0f * dif * t;//(*mobject)[i].GetColor().b;
                }
            }
            i = 0;
            xcoord += (xstep * 4);
            index += 4;
        }
        ycoord += ystep;
    }
}

#endif